#include <iostream>
#include <iomanip>
#include <math.h>
#include <graphics.h>

using namespace std;

float dt=1e-2;
float pi=3.1415926, k=2e-1, diss=0.2;

class vector{
  private:
    float coords[3];
  public:
    vector(float x=0,float y=0,float z=0){
	  coords[0] = x;
	  coords[1] = y;
	  coords[2] = z;
	}
	vector(vector & obj){
	  coords[0] = obj.coords[0];
	  coords[1] = obj.coords[1];
	  coords[2] = obj.coords[2];
	}
	void set(float x=0,float y=0,float z=0){
	  coords[0] = x;
	  coords[1] = y;
	  coords[2] = z;
	}
	float getx(void){ return coords[0]; }	
	float gety(void){ return coords[1]; }
	float getz(void){ return coords[2]; }
	void stretch(float koeff){
	  coords[0] *= koeff;
	  coords[1] *= koeff;
	  coords[2] *= koeff;
	}
	float length(void){
	  float len = pow(coords[0],2);
	  len += pow(coords[1],2);
	  len += pow(coords[2],2);
	  return len;
	}
	friend vector& summ(vector& obj1, vector& obj2){
	  static vector tempplus;
	  tempplus.coords[0]=obj1.coords[0]+obj2.coords[0];	  
	  tempplus.coords[1]=obj1.coords[1]+obj2.coords[1];
	  tempplus.coords[2]=obj1.coords[2]+obj2.coords[2];
	  return tempplus;
	}
	friend vector& operator - (vector& obj1){
	  static vector tempplus;
	  tempplus.coords[0]=-obj1.coords[0];	  
	  tempplus.coords[1]=-obj1.coords[1];
	  tempplus.coords[2]=-obj1.coords[2];
	  return tempplus;
	}
	friend vector& operator * (vector& obj1, float k){
	  static vector tempmult;
      tempmult.coords[0]=obj1.coords[0]*k;   
      tempmult.coords[1]=obj1.coords[1]*k;
      tempmult.coords[2]=obj1.coords[2]*k;
      return tempmult;
	}
	friend float scalar(vector& obj1, vector& obj2){
	  float result=0;
	  for(int i=0;i<=2;i++) result+= (obj1.coords[i]*obj2.coords[i]);
	  return result;
	}
	void print(void){
	  cout << "Vector output\n" << fixed<<setprecision(2);
	  cout << coords[0] << endl;
	  cout << coords[1] << endl;
	  cout << coords[2] << endl;
	}
	
};


class point{
  private:
    float phi, theta;
	float v_phi, v_theta;
	float a_phi, a_theta;
  public:    	
	vector r, temp, F;
	vector e_phi, e_theta;
	
    void set(float p, float t){
	  phi = p;
	  theta = t;
	  v_phi = v_theta = a_phi = a_theta = 0;
	  r.set(sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta));
    }	
	void print(void){
	  cout << phi << endl;
	  cout << theta << endl << endl;
	}
	void force(point& p){        //?????????????? c ???????? p
	  temp=summ(r,-(p.r));       //????????????? ??????????
	  temp.stretch(k/pow(temp.length(),3)); //??????????? ???? ?? ??????? ??????? p
	  F=summ(F,temp);            // ???????? ?????????? ????
    }
    void move(void){
	  e_theta.set(cos(theta)*cos(phi), cos(theta)*sin(phi), -sin(theta));
	  e_phi.set(-sin(phi), cos(phi), 0);    //???? ????????????? ???????
	  
	  a_theta=scalar(F,e_theta);            //???????????? ????????? ??? ?????????
	  a_phi=(scalar(F,e_phi) - v_theta*cos(phi)*v_phi)/sin(theta);
	
	  v_phi+=a_phi*dt;      //?????????? ?????????
	  v_theta+=a_theta*dt;
	  
	  temp.set(sin(theta)*v_phi,v_theta,0);
	  temp.stretch(1-diss*dt);
	  v_phi=temp.getx()/sin(theta);
	  v_theta=temp.gety();
	  
	  
	  phi+=v_phi*dt;        //?????????? ?????????
	  theta+=v_theta*dt;
	  r.set(sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta));
	}
	
};


int main(){
  
  int N=6, frames=25;
  int i,j, ticks;
  int front_x=300, front_y=125, top_x=300, top_y=425;
  
  point ball[N];
    ball[1].set(0,pi/2);       //????????? ????????? ??????
	ball[0].set(pi/4,pi/4);
	ball[2].set(pi/2,-pi/3);
	ball[3].set(-pi/8,pi/2);
	ball[4].set(pi/6,-pi/4);
	ball[5].set(4*pi/3,-pi/7);
	
	ticks = 0;
        
 
  initwindow(600, 550,"workspace");
  setcolor(14);
  circle(front_x,front_y,100);  
  circle(top_x,top_y,100);
  getch();
 
	while(true){
	
	  ticks++;
	
	  for(i=0;i<N;i++){                //???? ???????????????
	    ball[i].F.set(0,0,0);
	    for(j=0;j<N;j++){
		  if(j==i) continue;
          ball[i].force(ball[j]);
		}  
	  }
	  
	  
	  for(i=0;i<N;i++) ball[i].move();  //???? ?????????	  
	  
	  if(ticks > frames){
	    ticks = 0;
		setcolor(0); 
        setfillstyle(1,0);
		fillellipse(front_x,front_y,110,110);
		fillellipse(top_x,top_y,110,110);
		setcolor(14);
        circle(front_x,front_y,100);  
        circle(top_x,top_y,100);
		
		for(i=0;i<N;i++){
		  setcolor(i+1);
		  setfillstyle(1,i+1);
		  fillellipse(int(ball[i].r.getx()*100)+top_x,-int(ball[i].r.gety()*100)+top_y,5,5);
		  fillellipse(int(ball[i].r.getx()*100)+front_x,-int(ball[i].r.getz()*100)+front_y,5,5);
		  delay(30);
		}
		
	  }
	  
	}
	
 
  
  return 0;
}
